8  scTCRseq: Plot clone frequency over time

8.1 Set up workspace

# Libraries
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp

Attaching package: 'SeuratObject'

The following objects are masked from 'package:base':

    intersect, t

library(dplyr)
library(scRepertoire)
library(patchwork)
library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:patchwork':

    align_plots

The following object is masked from 'package:lubridate':

    stamp

library(grid)
library(forcats)

8.2 Load P101 and P103 TILs

p101_tils <- read.csv("sctcr_scRep_p101_til_freq_Part2.csv")
p103_tils <- read.csv("sctcr_scRep_p103_til_freq_Part2.csv")

8.3 Format for frequency plots

p101_tils_long <- p101_tils %>%
  pivot_longer(cols=c("P101_Tumor_W00_freq", "P101_Tumor_W12_freq", "P101_Tumor_W20_freq"), values_to = "Freq", names_to = "Timepoint") %>%
  mutate(patient = "P101") %>%
  select(-P101_Tumor_W00_count, -P101_Tumor_W12_count, -P101_Tumor_W20_count)
p103_tils_long <- p103_tils %>%
  pivot_longer(cols=c("P103_Tumor_W00_freq", "P103_Tumor_W12_freq", "P103_Tumor_W20_freq"), values_to = "Freq", names_to = "Timepoint") %>%
  mutate(patient = "P103") %>%
  select(-P103_Tumor_W00_count, -P103_Tumor_W12_count, -P103_Tumor_W20_count)

tils_long <- rbind(p101_tils_long, p103_tils_long)

tils_long <- tils_long %>%
  mutate(sctcr_category = factor(sctcr_category, levels = c("Existing", "Post-Nivolumab", "Post-Vaccine")),
         sctcr_category2 = factor(sctcr_category2, levels = c("Existing", "Post-Nivolumab", "Post-Vaccine", "Reactive")),
         sctcr_category3 = factor(sctcr_category3, levels = c("Existing", "Post-Nivolumab", "Post-Vaccine", "Tested, non-reactive", "Reactive")),
         vjaa = forcats::fct_reorder(as.factor(vjaa), rank(sctcr_category)),
         Timepoint = str_split_i(Timepoint, "_", 3))
  # mutate(Freq = case_when(Freq == 0 ~ NA,
  #                         Freq > 0 ~ Freq))

8.4 Summarize number of clones per category

tils_long %>%
  distinct(patient, vjaa, sctcr_category) %>%
  dplyr::count(patient, sctcr_category) %>%
  pivot_wider(names_from = c("patient"), values_from = "n")

# A tibble: 3 × 3
  sctcr_category  P101  P103
  <fct>          <int> <int>
1 Existing         300   516
2 Post-Nivolumab   810   890
3 Post-Vaccine     246   612

tils_long %>%
  distinct(patient, vjaa, sctcr_category2) %>%
  dplyr::count(patient, sctcr_category2) %>%
  pivot_wider(names_from = c("patient"), values_from = "n")

# A tibble: 4 × 3
  sctcr_category2  P101  P103
  <fct>           <int> <int>
1 Existing          300   516
2 Post-Nivolumab    810   890
3 Post-Vaccine      243   609
4 Reactive            3     3

tils_long %>%
  distinct(patient, vjaa, sctcr_category3) %>%
  dplyr::count(patient, sctcr_category3) %>%
  pivot_wider(names_from = c("patient"), values_from = "n")

# A tibble: 5 × 3
  sctcr_category3       P101  P103
  <fct>                <int> <int>
1 Existing               297   513
2 Post-Nivolumab         807   884
3 Post-Vaccine           233   607
4 Tested, non-reactive    16    11
5 Reactive                 3     3

8.5 Create a pie chart of number of clones per category

pie <- tils_long %>%
  distinct(patient, vjaa, sctcr_category) %>%
  dplyr::count(patient, sctcr_category) %>%
  ggplot(aes(x="", y=n, fill=sctcr_category)) +
  facet_wrap(~patient, ncol = 2, scales = "free") +
  geom_bar(stat="identity", width=1, color = "black", size = 1) +
  coord_polar("y", start=0) +
  theme_void() +
  labs(fill = "TCR type") +
  scale_fill_manual(values = c("grey80", "darkgoldenrod2", "deepskyblue3")) +
  theme(strip.text = element_text(size = 20)) +
  geom_text(aes(label = n),
            position = position_stack(vjust = 0.5), 
            size=4, 
            color = "black", 
            fontface = "bold")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

pie

8.6 Frequency plot for Fig 4G

# Line graphs
p101_line <- tils_long %>%
  filter(patient == "P101") %>%
  ggplot(aes(x = Timepoint, y = Freq, color = sctcr_category, group = vjaa)) +
  geom_line(position = position_jitter(width = 0, height = 0.1, seed = 123), size = 0.35) +
            # , aes(size = sctcr_category2)) +
  facet_wrap(~patient, scale = "free_x") +
  theme_classic() +
  scale_color_manual(values = c("grey80", "darkgoldenrod2", "deepskyblue3")) +
  theme(axis.text.x = element_text(face ="bold", size = 12), 
        strip.text = element_text(face="bold", size=12),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)) +
  scale_y_log10(labels = scales::comma) +
  labs(color = "TCR type") +
  ylab("log10 Clonotype Frequency (jittered)") +
  annotation_logticks(sides = "l", scaled = TRUE, linewidth = 0.2)
  # coord_cartesian(ylim = c(0.0003, 0.1), expand = TRUE)

p103_line <- tils_long %>%
  filter(patient == "P103") %>%
  ggplot(aes(x = Timepoint, y = Freq, color = sctcr_category, group = vjaa)) +
  geom_line(position = position_jitter(width = 0, height = 0.1, seed = 123), size = 0.35) +
            # , aes(size = sctcr_category2)) +
  facet_wrap(~patient, scale = "free_x") +
  theme_classic() +
  scale_color_manual(values = c("grey80", "darkgoldenrod2", "deepskyblue3")) +
  labs(color = "TCR type") +
  ylab("log10 Clonotype Frequency (jittered)") +
  theme(axis.text.x = element_text(face ="bold", size = 12), 
        strip.text = element_text(face="bold", size=12),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)) +
  scale_y_log10(labels = scales::comma, , breaks = NULL)
  # coord_cartesian(ylim = c(0.0003, 0.1), expand = TRUE)

(p101_line + p103_line) +
  plot_layout(guides = "collect", 
              axis_titles = "collect") & theme(legend.position = "right", 
                                               legend.title=element_text(size=15),
                                               legend.text = element_text(size = 15))

Warning in scale_y_log10(labels = scales::comma): log-10 transformation
introduced infinite values.

Warning in scale_y_log10(labels = scales::comma, , breaks = NULL): log-10
transformation introduced infinite values.

(p101_line + p103_line) +
  plot_layout(guides = "collect", 
              axis_titles = "collect") & theme(legend.position = "right", 
                                               legend.title=element_text(size=15),
                                               legend.text = element_text(size = 15))

Warning in scale_y_log10(labels = scales::comma): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.

8.7 Get session info

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] cowplot_1.1.3      patchwork_1.3.0    scRepertoire_2.0.0 Seurat_5.1.0      
 [5] SeuratObject_5.0.2 sp_2.2-0           lubridate_1.9.4    forcats_1.0.0     
 [9] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.4        readr_2.1.5       
[13] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] cubature_2.1.1              RcppAnnoy_0.0.22           
  [3] splines_4.3.2               later_1.4.1                
  [5] bitops_1.0-9                polyclip_1.10-7            
  [7] fastDummies_1.7.5           lifecycle_1.0.4            
  [9] globals_0.16.3              lattice_0.22-7             
 [11] MASS_7.3-60.0.1             magrittr_2.0.3             
 [13] plotly_4.10.4               rmarkdown_2.29             
 [15] httpuv_1.6.15               sctransform_0.4.1          
 [17] spam_2.11-1                 spatstat.sparse_3.1-0      
 [19] reticulate_1.42.0           pbapply_1.7-2              
 [21] RColorBrewer_1.1-3          abind_1.4-8                
 [23] zlibbioc_1.48.2             Rtsne_0.17                 
 [25] GenomicRanges_1.54.1        ggraph_2.2.1               
 [27] BiocGenerics_0.48.1         RCurl_1.98-1.17            
 [29] tweenr_2.0.3                evmix_2.12                 
 [31] GenomeInfoDbData_1.2.11     IRanges_2.36.0             
 [33] S4Vectors_0.40.2            ggrepel_0.9.5              
 [35] irlba_2.3.5.1               listenv_0.9.1              
 [37] spatstat.utils_3.1-0        iNEXT_3.0.1                
 [39] MatrixModels_0.5-3          goftest_1.2-3              
 [41] RSpectra_0.16-2             spatstat.random_3.3-1      
 [43] fitdistrplus_1.2-2          parallelly_1.41.0          
 [45] leiden_0.4.3.1              codetools_0.2-20           
 [47] DelayedArray_0.28.0         ggforce_0.4.2              
 [49] tidyselect_1.2.1            farver_2.1.2               
 [51] viridis_0.6.5               matrixStats_1.5.0          
 [53] stats4_4.3.2                spatstat.explore_3.3-2     
 [55] jsonlite_1.8.9              tidygraph_1.3.1            
 [57] progressr_0.15.1            ggridges_0.5.6             
 [59] ggalluvial_0.12.5           survival_3.8-3             
 [61] tools_4.3.2                 stringdist_0.9.12          
 [63] ica_1.0-3                   Rcpp_1.0.14                
 [65] glue_1.8.0                  gridExtra_2.3              
 [67] SparseArray_1.2.4           xfun_0.50                  
 [69] MatrixGenerics_1.14.0       GenomeInfoDb_1.38.8        
 [71] withr_3.0.2                 fastmap_1.2.0              
 [73] SparseM_1.84-2              digest_0.6.37              
 [75] timechange_0.3.0            R6_2.6.1                   
 [77] mime_0.13                   colorspace_2.1-1           
 [79] scattermore_1.2             tensor_1.5                 
 [81] spatstat.data_3.1-2         utf8_1.2.4                 
 [83] generics_0.1.3              data.table_1.15.4          
 [85] graphlayouts_1.1.1          httr_1.4.7                 
 [87] htmlwidgets_1.6.4           S4Arrays_1.2.1             
 [89] uwot_0.2.3                  pkgconfig_2.0.3            
 [91] gtable_0.3.6                lmtest_0.9-40              
 [93] SingleCellExperiment_1.24.0 XVector_0.42.0             
 [95] htmltools_0.5.8.1           dotCall64_1.2              
 [97] scales_1.3.0                Biobase_2.62.0             
 [99] png_0.1-8                   spatstat.univar_3.0-0      
[101] ggdendro_0.2.0              knitr_1.49                 
[103] rstudioapi_0.17.1           rjson_0.2.23               
[105] tzdb_0.5.0                  reshape2_1.4.4             
[107] nlme_3.1-168                zoo_1.8-13                 
[109] cachem_1.1.0                KernSmooth_2.23-26         
[111] parallel_4.3.2              miniUI_0.1.1.1             
[113] pillar_1.10.1               vctrs_0.6.5                
[115] RANN_2.6.2                  VGAM_1.1-13                
[117] promises_1.3.2              xtable_1.8-4               
[119] cluster_2.1.8.1             evaluate_1.0.1             
[121] truncdist_1.0-2             cli_3.6.3                  
[123] compiler_4.3.2              rlang_1.1.5                
[125] crayon_1.5.3                future.apply_1.11.3        
[127] labeling_0.4.3              plyr_1.8.9                 
[129] stringi_1.8.4               viridisLite_0.4.2          
[131] deldir_2.0-4                munsell_0.5.1              
[133] gsl_2.1-8                   lazyeval_0.2.2             
[135] spatstat.geom_3.3-2         quantreg_6.1               
[137] Matrix_1.6-5                RcppHNSW_0.6.0             
[139] hms_1.1.3                   future_1.34.0              
[141] shiny_1.9.1                 SummarizedExperiment_1.32.0
[143] evd_2.3-7.1                 ROCR_1.0-11                
[145] igraph_2.0.3                memoise_2.0.1